function udt = h2eom(t, u, params)
%   choice of Pi: 
%   Pi_xx =  2 P1 (1+UX^2) +  2 UX UY P2
%   Pi_yy = -2 P1 (1+UY^2) +  2 UX UY P2
%   Pi_xy =  P2 ( 2 + UX^2 + UY^2 )

    A = params.A;
    Nx = params.Nx;
    Ny = params.Ny;
    Nt = Nx.*Ny;

    %u = gpuArray(uu);
    


    % unpack variables 
    T  = reshape(u(0*Nt+1 : 1*Nt), Nx, Ny);
    ux = reshape(u(1*Nt+1 : 2*Nt), Nx, Ny);
    uy = reshape(u(2*Nt+1 : 3*Nt), Nx, Ny);
    P1 = reshape(u(3*Nt+1 : 4*Nt), Nx, Ny);
    P2 = reshape(u(4*Nt+1 : 5*Nt), Nx, Ny);

    Tdx  = dx(T);
    Tdy  = dy(T);
    uxdx = dx(ux);
    uxdy = dy(ux);
    uydx = dx(uy);
    uydy = dy(uy);
    P1dx = dx(P1);
    P1dy = dy(P1);
    P2dx = dx(P2);
    P2dy = dy(P2);


    
    Tdt  = (1/6).*T.^(-2).*(3.*T.^2.*(Tdx.*ux+Tdy.*uy)+4.*((-1).*P1dx.*ux+(-1).*P2dy.*ux+P1dy.*uy+(-1).*P2dx.*uy+P2.*(uxdy+uydx))+4.*P1.*(uxdx+(-1).*uydy)+3.*T.^3.*(uxdx+uydy));

    uxdt = (16.*(P1.^2+P2.^2)+(-9).*T.^6).^(-1).*(16.*P2.*((-1).*P1dy+P2dx)+(4.*P1+(-3).*T.^3).*(4.*(P1dx+P2dy)+3.*T.^2.*Tdx)+12.*P2.*T.^2.*Tdy)+(1/6).*(1+A).^(-1).*T.^(-3).*((-16).*(P1.^2+P2.^2)+9.*T.^6).^(-1).*(128.*(1+A).*P2.^3.*ux.*(uxdy+uydx)+9.*T.^9.*(4.*uxdy.*uy+6.*A.*uxdy.*uy+(-2).*uy.*uydx+(1+3.*A).*ux.*(uxdx+(-1).*uydy))+128.*(1+A).*P1.^3.*ux.*(uxdx+(-1).*uydy)+48.*P2.^2.*T.^3.*(2.*T.*ux+((-1)+A).*uxdy.*uy+(-1).*(1+3.*A).*uy.*uydx+(1+A).*ux.*(uxdx+uydy))+(-12).*P2.*T.^6.*(ux.*(uxdy+(-3).*A.*uxdy+(7+15.*A).*uydx)+2.*uy.*(3.*T+uxdx+2.*uydy+3.*A.*uydy))+32.*P1.^2.*(3.*T.^4.*ux+4.*(1+A).*P2.*ux.*(uxdy+uydx)+3.*T.^3.*(A.*uy.*(uxdy+(-1).*uydx)+(1+A).*ux.*uydy))+(-4).*P1.*(18.*T.^7.*ux+12.*(1+A).*P2.*T.^3.*(ux.*(uxdy+uydx)+uy.*(uxdx+(-1).*uydy))+(-32).*(1+A).*P2.^2.*ux.*(uxdx+(-1).*uydy)+3.*T.^6.*(2.*(1+3.*A).*uy.*(2.*uxdy+(-1).*uydx)+ux.*(7.*uxdx+9.*A.*uxdx+(-1).*uydy+(-3).*A.*uydy))));

    uydt = (16.*(P1.^2+P2.^2)+(-9).*T.^6).^(-1).*(4.*(4.*P2.*(P1dx+P2dy)+(P1dy+(-1).*P2dx).*(4.*P1+3.*T.^3)+3.*P2.*T.^2.*Tdx)+(-3).*T.^2.*(4.*P1+3.*T.^3).*Tdy)+(1/6).*(1+A).^(-1).*T.^(-3).*((-16).*(P1.^2+P2.^2)+9.*T.^6).^(-1).*(128.*(1+A).*P2.^3.*uy.*(uxdy+uydx)+16.*P2.^2.*(6.*T.^4.*uy+8.*(1+A).*P1.*uy.*(uxdx+(-1).*uydy)+3.*T.^3.*((-1).*ux.*(uxdy+3.*A.*uxdy+uydx+(-1).*A.*uydx)+(1+A).*uy.*(uxdx+uydy)))+4.*P2.*((-18).*T.^7.*ux+32.*(1+A).*P1.^2.*uy.*(uxdy+uydx)+12.*(1+A).*P1.*T.^3.*(uy.*(uxdy+uydx)+ux.*((-1).*uxdx+uydy))+(-3).*T.^6.*(uy.*(7.*uxdy+15.*A.*uxdy+uydx+(-3).*A.*uydx)+2.*ux.*(2.*uxdx+3.*A.*uxdx+uydy)))+(4.*P1+3.*T.^3).*(24.*P1.*T.^4.*uy+32.*P1.^2.*uy.*(uxdx+(-1).*uydy)+24.*P1.*T.^3.*uy.*uydy+3.*T.^6.*((-2).*ux.*uxdy+(-1).*uxdx.*uy+4.*ux.*uydx+uy.*uydy)+A.*(32.*P1.^2.*uy.*(uxdx+(-1).*uydy)+24.*P1.*T.^3.*((-1).*ux.*uxdy+ux.*uydx+uy.*uydy)+9.*T.^6.*((-1).*uxdx.*uy+2.*ux.*uydx+uy.*uydy))));

    P1dt = (1/12).*(12.*P1dx.*ux+12.*P1dy.*uy+(1+A).^(-1).*T.^(-3).*((-16).*(P1.^2+P2.^2)+9.*T.^6).^(-1).*(4.*(P1dx+P2dy)+3.*T.^2.*Tdx).*((64.*(1+A).*P1.*(P1.^2+P2.^2)+24.*((1+A).*P1.^2+(-2).*P2.^2).*T.^3+(-6).*(7+9.*A).*P1.*T.^6+(-9).*T.^9).*ux+24.*P2.*T.^3.*((3+A).*P1+(-2).*T.^3).*uy)+(-1).*(1+A).^(-1).*T.^(-3).*((-16).*(P1.^2+P2.^2)+9.*T.^6).^(-1).*(4.*P1dy+(-4).*P2dx+(-3).*T.^2.*Tdy).*(24.*P2.*T.^3.*((3+A).*P1+2.*T.^3).*ux+(64.*(1+A).*P1.*(P1.^2+P2.^2)+(-24).*((1+A).*P1.^2+(-2).*P2.^2).*T.^3+(-6).*(7+9.*A).*P1.*T.^6+9.*T.^9).*uy)+(1+A).^(-1).*T.^(-3).*(3.*T.^3.*(4.*A.*P2.*((-1).*uxdy+uydx)+T.^3.*(uxdx+(-1).*uydy))+16.*(1+A).*P1.^2.*(uxdx+(-1).*uydy)+2.*P1.*(6.*T.^4+8.*(1+A).*P2.*(uxdy+uydx)+9.*(1+A).*T.^3.*(uxdx+uydy))));

    P2dt = (1/12).*(12.*P2dx.*ux+12.*P2dy.*uy+(1+A).^(-1).*T.^(-3).*((-16).*(P1.^2+P2.^2)+9.*T.^6).^(-1).*(4.*(P1dx+P2dy)+3.*T.^2.*Tdx).*(2.*P2.*(32.*(1+A).*(P1.^2+P2.^2)+12.*(3+A).*P1.*T.^3+(-3).*(7+9.*A).*T.^6).*ux+3.*T.^3.*((-16).*P1.^2+8.*(1+A).*P2.^2+16.*P1.*T.^3+(-3).*T.^6).*uy)+(1+A).^(-1).*T.^(-3).*((-16).*(P1.^2+P2.^2)+9.*T.^6).^(-1).*(4.*P1dy+(-4).*P2dx+(-3).*T.^2.*Tdy).*(3.*T.^3.*((-8).*(1+A).*P2.^2+(4.*P1+T.^3).*(4.*P1+3.*T.^3)).*ux+2.*P2.*((-32).*(1+A).*(P1.^2+P2.^2)+12.*(3+A).*P1.*T.^3+3.*(7+9.*A).*T.^6).*uy)+(1+A).^(-1).*T.^(-3).*(16.*(1+A).*P2.^2.*(uxdy+uydx)+3.*T.^3.*(4.*A.*P1.*(uxdy+(-1).*uydx)+T.^3.*(uxdy+uydx))+2.*P2.*(6.*T.^4+8.*(1+A).*P1.*(uxdx+(-1).*uydy)+9.*(1+A).*T.^3.*(uxdx+uydy))));



    udt = zeros(Nx.*Ny.*5, 1);

    udt(0*Nt+1 : 1*Nt) = Tdt(:);
    udt(1*Nt+1 : 2*Nt) = uxdt(:);
    udt(2*Nt+1 : 3*Nt) = uydt(:);
    udt(3*Nt+1 : 4*Nt) = P1dt(:);
    udt(4*Nt+1 : 5*Nt) = P2dt(:);

    % udt(0*Nt+1 : 1*Nt) = gather(Tdt(:));
    % udt(1*Nt+1 : 2*Nt) = gather(uxdt(:));
    % udt(2*Nt+1 : 3*Nt) = gather(uydt(:));
    % udt(3*Nt+1 : 4*Nt) = gather(P1dt(:));
    % udt(4*Nt+1 : 5*Nt) = gather(P2dt(:));


end
